Intriguing Heat Conduction of a Polymer Chain 
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We study heat conduction in a one- dimensional (ID) chain of particles with longitudinal as well 
as transverse motions. The particles are connected by two-dimensional harmonic springs together 
with bending angle interactions. Using equilibrium and nonequilibrium molecular dynamics, three 
types of thermal conducting behaviors are found: a logarithmic divergence with system sizes for large 
transverse coupling, 1/3 power- law at intermediate coupling, and 2/5 power-law at low temperatures 
and weak coupling. The results are consistent with a simple mode-coupling analysis of the same 
model. The 1/3 power- law divergence should be a generic feature for models with transverse motions. 
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To understand the microscopic dynamical mechanism 
of heat conduction is one of the long standing tasks in 
nonequilibrium statistical mechanics. This problem has 
attracted increasing attention in recent years [1-13]. The 
main effort has been focused on the necessary and suf- 
ficient conditions of the Fourier law of heat conduction. 
With strong numerical support, it is argued that chaos 
(or exponential instability) is the necessary condition . 
However, recent results show that even linear instability, 
such as that found in generic polygonal billiards, is suffi- 
cient for a normal diffusion and energy transport obeying 
the Fourier law 0. 

In many systems studied so far, the heat conduction 
violates the Fourier law, namely, the thermal conduc- 
tivity diverges with system size N as iV", with a > 0. 
It has been proved that for ID system, the momentum 
conservation leads to an anomalous heat conduction 
However, its specific form of divergence with system size 
is still of considerable controversy 4, 5, 6, 7, 8J. Based on 
a renormalization group analysis for a ID hydrodynamic 
fluid model, it is argued that in a generic momentum 
conserving system, the thermal conductivity should be 
K oc N^^^Ai\. Unfortunately, most existing numerical re- 
sults do not agree with this prediction because all ID 
lattice models considered so far have no transverse de- 
gree of freedom that is required in the analysis. However, 
a mode-coupling theory analysis for the ID Fermi-Pasta- 
Ulam (FPU) model gives a divergent exponent 2/5 0,IH, 
which is supported by the numerics from different groups 
0. 

It seems that a universal exponent does not exist. Most 
recently, it is found that a divergent thermal conductivity 
is connected with a supcrdiffusion k oc iV^"^/'', where 
f3 is the exponent of the diffusion (Aa;^ t^, < /? < 2). 
The value of f3 changes from model to model. This is 
justified by all billiard gas models studied. 

On the other hand, the understanding of heat conduc- 



tion mechanism will allow us to control and manipulate 
heat current, and eventually to design novel thermal de- 
vices with certain function [lo| . To this end, more re- 
alistic physical models are necessary. Among many oth- 
ers, nanotubes and polymer chains are most promising. 
Recent molecular dynamics (MD) study of Carbon nan- 
otubes with realistic interaction potential suggested a di- 
vergent thermal conductivity for narrow diameter tubes 
[m . The strict ID models may not be applicable to nan- 
otubes. The added transverse motion and the flexibility 
of the tube at long length scales will certainly scatter 
phonons, and thus should have a profound effect on ther- 
mal transport. 

In this Letter, we consider a polymer chain of N point 
particles with mass m on a ID lattice. The lattice fixes 
the connectivity topology such that only the neighboring 
particles interact. The Hamiltonian is given by 
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where the position vector r = {x,y) and momentum vec- 
tor p — {px,Py) are two-dimensional; a is lattice con- 
stant. If the system is restricted to yi = 0, it is essen- 
tially a ID gas with harmonic interaction. The coupling 
Kr is the spring constant; signifies bending or flexi- 
bility of the chain, while (pi is the bond angle formed with 
two neighbor sites, cos 4>i = — nj_i • n^, and unit vector 
Ui = Ari/|Ar,|, Ar., = r,+i - r,. 

We determine the heat current in a temperature gra- 
dient by nonequilibrium MD. The system is set up with 
fixed left most and right most boundary. The average dis- 
tance between particles is set to a, the zero-temperature 
equilibrium distance. A group of four particles at the 
two ends are subject to heat baths at temperature 
and Th, respectively. This is realized by Nose- Hoover 
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rent multiplied by N , jN = {Th — in log-log 

scale (Fig. ^a)) and linear-log scale (Fig. H^b)). It is 
clearly shown that three types of behaviors of the ther- 
mal conductivity k are observed, the logarithmic diver- 
gence, log A^, power-law k oc N'^ with a = 1/3, as well as 
2/5, depending on the model parameters. Log-log plot 
shows linear behavior for data set E, F, H, and J. At 
the parameters of set E, excellent power-law dependence 
is found, with an exponent of a = 0.334 ± 0.003 (using 
an error weighted least-squares fit for N > 128). Set F 
is also in good agreement with a slope of 1/3. On the 
other hand, for set H and J, we have exponent a con- 
sistent with 0.4. Set B is consistent with logarithmic 
divergence, k oc logA^ (see Fig. ^b)). The model has 
two key parameters, the temperature T, and the trans- 
verse coupling K^. We should mention that wide range 
of parameters is scanned, and surprisingly, only the three 
scalings are found so far in this model. 

To understand the simulation results, we consider 
a simple mode-coupling theory for the present model. 
The equations of motion in terms of normal-mode co- 
ordinates, qI = v^Ef=o'(^. " ja)e'2.,'c/iv^ q± ^ 
-1 



y^^^"^^ j/j-e*^^-''^/^, for small oscillation near zero- 
temperature equilibrium position, keeping only leading 
nonlinearity, are 



FIG. 1: jN vs A'' on double-logarithmic (a) and linear-log 
plot (b). The parameters {K^,Tl,Th) of the model are, set 
B: (1, 0.2, 0.4), set E: (0.3, 0.3, 0.5), set F: (1, 5, 7), set H: 
(0, 0.3, 0.5), set J: (0.05, 0.1, 0.2). AU of them have Kr = 1, 
mass m = 1, lattice spacing a — 2. The straight lines on F 
and E have slope 1/3, while the slope on H and J is 2/5. 



^ = -i-b'Ql + E 4,,''Qi'Qi'. (3a) 

k' + k" = k 

- = -i^irQi+ J2 4',k"QlQb', (3b) 

k'+k"=k 



thermostats. The rest of the particles follow the equa- 
tions of motion using a velocity Verlet algorithm. We 
use small time step sizes h = 0.0005 to 0.0010. Typical 
MD steps are 10® to lO^"^. 

We use the following expression for local heat current 
per particle: 

mji = - Ari((p,j + Pj+i) • G(i)) 

- Ar,_i((p, +p,_i) • G{i - 1)) 
-t-Ar,_i(p, •H(^-2,^-l,^-l)) 
-t-Ari(p, • H(i + l,i-|-l,i)) -|-p,/ii, (2) 

where G{i) = \Kr{\l^ri\ - a)ni, ll{i,j,k) = K^{ni + 
rifc cos(/)j)/|Arfe|, and the local energy per particle hi — 
lKr{i\Ar,_i\-a)^ + i\Ar,\-a)^)+K^cos(l),+pf/{2m). 
This is derived from J = '^d[rihi) / dt, by regrouping 
some of the terms using translational invariance. It sat- 
isfies the continuity equation in the long- wave limit. 

We present our main numerical results in Fig. ^ The 
data are obtained using 20 IGHz-Pentium PCs over six 
months of CPU times. We plot the average heat cur- 



where the bare dispersion relations are given by a;[: = 
2^/^|sin^|, and ui = 'i\[^sm^^. The expres- 



Trk 

N 



sions for c):' are complicated, but can be simplified in 

the long-wave limit, as c^^,, — 2c| j,, oc kk'{k + k'). In- 
stead of the integer fc, we can also index the mode by its 
corresponding lattice momentum p — 2TTk/{aN). 

A central quantity in the mode-coupling the- 
ory is the normalized correlation function, gp{t) = 
{Qp(t)Qp{0))/{\QpiQ)\^)- The Fourier-Laplace transform 



of the correlation function, 
given by [ll[li. 



9[z\ 
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■ dt, is 
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The constants cy and ca_ are effective or renormal- 
ized sound velocities for the longitudinal and transverse 
modes. They are defined, e.g., by (c||p)^(|(3|,p) — ksT, 
as p — > 0. The damping functions (memory kernel) are 
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given by the coupled equations, 

= 5 j\dp{9^it))\ (5a) 

= ^ j\dpgl{t)g^{t). (5b) 

Eqs. 101) and lO form a closed system of nonlinear in- 
tegral equations. This is a straightforward generaliza- 
tion of the strict ID result The above equations are 
derived under a number of simplification assumptions, 
such as long-wave approximation, mean-field type prod- 
uct approximation for the correlation functions, replac- 
ing random-force correlation with true force correlation. 
Some of them can be removed but more complicated 
equations will result. 

In Fourier space, for large z the solution is found from 
integration by part, as j/ll'-'-[z] = j_/(iza) + 0{z^'^). 
The long-wave asymptotic decay of each mode is con- 
trolled by the small z behavior of the function 
We define <5|| and 5±^ by oc z^^w-^. The dispersion 

relation is then given by the location of the poles in the 
correlation function g[z\. The imaginary part of the fre- 
quency gives damping, by 7p oc p'^v[z\\^_^^^ oc p^^^ . We 
note that three types of behaviors can be derived from 
the above set of equations. If K\\ « K_\_ and c\\ ~ c_l, 
the two equations reduce to that of strict ID model, we 
thus expect the result of Lepri |^, i.e., 5\\ = 5^ = 1/3. 
One the other hand, it can be shown rigorously that in 
the limit of small and small c_l , we have = and 

= 1/2. Formally, when a ^ 0, the equation possesses 
the scaling solution of the form v\\z] = A~"'^i^[2]; this im- 
plies v[z\ oc 1/z. These analytic results are supported by 
numerical solutions of the coupled equations, shown in 
Fig. 121 They are solved by a brute force numerical in- 
tegration in Fourier space. Details of the mode-coupling 
calculation will be presented elsewhere. 

In Fig. 13 at parameter set I, we observe very good 
asymptotic behavior of v^^[z\ oc z~^/^ and v'^[z\ oc const. 
This corresponds to the behavior of MD results for data 
set E and F in Fig. When C|| = c±_ but K\^\ ^ Kj^ (set 
II), we see crossover from = 1/3 to 1/2. The curve III 
may be related to the logarithmic divergence. We note 
that a meaningful, direct mapping from the simulation 
parameters to mode-coupling parameters is not possible, 
due to qualitative nature of the theory. 

The prediction of i5|| = 1/2 and (5j^ = is checked 
against an equilibrium MD simulation in a microcanoni- 
cal ensemble with periodic boundary condition. We com- 
pute the normal-mode correlation {Qp{t)Q*{0)) for each 
mode specified by the lattice momentum p = 27r/c/ (aN). 
The functions are oscillatory with an exponential decay, 
cos(a;i)e~^p*. The decay constants are obtained by fit- 
ting the maximum amplitude as a function of time. The 
results are presented in Fig. 13 Comparing with results 
from smaller and larger system sizes, effect of finite sizes 
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FIG. 2: Real part of u^lz] and I'^fz] vs z for parameters 
a = 1, = 1, and {K±, cy, c±) I: (0.3, 2, 1), II: (1.8, 1, 1), 
III: (2, 1, 0.5). 
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FIG. 3: The decay rate 7" (dots) and 7^ (triangles) vs k for 
the parameters set B, E, and J at equilibrium temperature 
T = [Tl + Th)/2. The number on the line indicates the slope 
of the straight line. The system size is A'^ = 1024. 

appears rather small at iV" = 1024. Excellent agreement 
with mode-coupling theory (7 oc k'^~^) is obtained for 
data set E. However, for data set B and J, the slopes are 
not consistent with either logarithmic divergence for n or 
2/5 law. This may be interpreted as that we are still not 
in the asymptotic regime. 

To connect the result of dampin g of the modes with 
thermal conductivity, it is noted [a, llal that each mode 
contributes to the thermal transport independently. Un- 
der the linear-response theory, the Green-Kubo formula 
relates the current-current correlation to the thermal con- 
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FIG. 4: The Green-Kubo integrand, {J{t)J{0))/N, vs time t. 
The parameters are the same as that in Fig. |3 The straight 
lines have slope —2/3 (on E) and —3/5 (on J). 

ductivity by 

The decay rates for J are assumed to be the same as that 
for Q, thus (J(i)J(O)) cx J2p^^'Pi~"fp'^)- The amphtude 
of the exponential decay is approximately independent 
of p. Converting the summation to integral, we have 
(J(t)J(O)) oc The thermal conductivity on a 

finite lattice is obtained by integrating over t to a time of 
0{N), Thus kn oc iVi-i/(2-5||) ^ TV". When J,, = 1/2, 
we have a ~ 1/3, and when 6\\ = 1/3, the exponent 
a = 2/5. 

The current-current correlation functions are presented 
in Fig. 01 for the parameters corresponding to data set B, 
E, and J in Fig.Q For data set J, a power-law dependence 
is in excellent agreement with the theoretical expectation 
t°'~^ with a = 2/5. For set E, the curve is a bit steeper 
than expected. This may be due to finite sizes. For set 
B, where logarithmic divergence is observed, we do not 
observe good power-law behavior in the correlation. 

We need to clarify the relationship between the three 
types of observed behaviors in the nonequilibrium MD 
results. From mode-coupling point of view, the 1/3-law 
is generic and robust, while a = 2/5 should be eventu- 
ally crossover to 1/3 at long length scales. However, such 
a crossover is not observed in MD data. The crossover 
effect can be argued for a more general setting. More gen- 
eral mode-coupling equations for a generic interaction po- 
tential consistent with the symmetry would have an ad- 
ditional term of ^ J^^ dp{gl{t))'^ for Eq. (|S3l; Eq. (|5bjl 
remains the same. Such a term can appear either from 



cubic or quartic nonlinearity in potential. Contribution 
from this extra term decays in time t faster than the per- 
pendicular component contribution. Thus, the asymp- 
totic result of J|| — 1/2 remains true. The same should 
be also true even if a chain is allowed to move in three 
dimensions. If the parameter is sufficiently large, we 
may see exponent close to 0.4. The logarithmic diver- 
gence is a bit difficult to interpret. It may require a more 
refined theory than the present naive mode-coupling the- 
ory. 

In summary, we have observed three different scalings 
in a ID polymer chain. When the transverse motion cou- 
ples with the longitudinal motion, the thermal conductiv- 
ity diverges with system size with a 1/3 power-law. This 
has been demonstrated with a very high precision nu- 
merical result and explained in terms of a mode-coupling 
theory. In the weak coupling regime, a 2/5 power-law is 
observed which is consistent with the results observed in 
the FPU model. In the case of strong transverse coupling, 
a logarithmic divergent law is observed. 
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